//
// Created by tyler on 4/13/17.
//

#include "quadrature/QuadratureRule.hpp"

YAFEL_NAMESPACE_OPEN

/*
void QuadratureRule::get_tetrahedron_quadrature(int porder)
{

    switch(porder) {
        case 1:
        default:
            nodes = {{1./4., 1./4., 1./4.}};
            weights = {1./6.};
            return;
    }


}
*/



void QuadratureRule::get_tetrahedron_quadrature(int porder)
{

    switch (porder) {
        case 1:
            nodes = {
                    {0.25, 0.25, 0.25},
            };
            weights = {0.16666666666666666};
            return;
        case 2:
            nodes = {
                    {0.138196601125011, 0.138196601125011, 0.585410196624969},
                    {0.138196601125011, 0.138196601125011, 0.138196601125011},
                    {0.138196601125011, 0.585410196624969, 0.138196601125011},
                    {0.585410196624969, 0.138196601125011, 0.138196601125011},
            };
            weights = {0.041666666666666664,
                       0.041666666666666664,
                       0.041666666666666664,
                       0.041666666666666664};
            return;
        case 3:
            nodes = {
                    {0.25,                0.25,                0.25},
                    {0.16666666666666666, 0.5,                 0.16666666666666666},
                    {0.16666666666666666, 0.16666666666666666, 0.5},
                    {0.16666666666666666, 0.16666666666666666, 0.16666666666666666},
                    {0.5,                 0.16666666666666666, 0.16666666666666666},
            };
            weights = {-0.13333333333333333,
                       0.075,
                       0.075,
                       0.075,
                       0.075};
            return;
        case 4:
            nodes = {
                    {0.25,              0.25,              0.25},
                    {0.071428571428571, 0.071428571428571, 0.785714285714286},
                    {0.071428571428571, 0.071428571428571, 0.071428571428571},
                    {0.071428571428571, 0.785714285714286, 0.071428571428571},
                    {0.785714285714286, 0.071428571428571, 0.071428571428571},
                    {0.399403576166799, 0.100596423833201, 0.100596423833201},
                    {0.100596423833201, 0.399403576166799, 0.399403576166799},
                    {0.100596423833201, 0.100596423833201, 0.399403576166799},
                    {0.399403576166799, 0.100596423833201, 0.399403576166799},
                    {0.100596423833201, 0.399403576166799, 0.100596423833201},
                    {0.399403576166799, 0.399403576166799, 0.100596423833201},
            };
            weights = {-0.013155555555555556,
                       0.007622222222222222,
                       0.007622222222222222,
                       0.007622222222222222,
                       0.007622222222222222,
                       0.024888888888888887,
                       0.024888888888888887,
                       0.024888888888888887,
                       0.024888888888888887,
                       0.024888888888888887,
                       0.024888888888888887};
            return;
        case 5:
            nodes = {
                    {0.25,                0.25,                0.25},
                    {0.3333333333333333,  0.3333333333333333,  0.3333333333333333},
                    {0.0,                 0.3333333333333333,  0.3333333333333333},
                    {0.3333333333333333,  0.0,                 0.3333333333333333},
                    {0.3333333333333333,  0.3333333333333333,  0.0},
                    {0.09090909090909091, 0.7272727272727273,  0.09090909090909091},
                    {0.09090909090909091, 0.09090909090909091, 0.7272727272727273},
                    {0.09090909090909091, 0.09090909090909091, 0.09090909090909091},
                    {0.7272727272727273,  0.09090909090909091, 0.09090909090909091},
                    {0.433449846426336,   0.066550153573664,   0.433449846426336},
                    {0.433449846426336,   0.433449846426336,   0.066550153573664},
                    {0.066550153573664,   0.066550153573664,   0.433449846426336},
                    {0.066550153573664,   0.433449846426336,   0.433449846426336},
                    {0.433449846426336,   0.066550153573664,   0.066550153573664},
                    {0.066550153573664,   0.433449846426336,   0.066550153573664},
            };
            weights = {0.030283678097089,
                       0.006026785714286,
                       0.006026785714286,
                       0.006026785714286,
                       0.006026785714286,
                       0.011645249086029,
                       0.011645249086029,
                       0.011645249086029,
                       0.011645249086029,
                       0.010949141561386,
                       0.010949141561386,
                       0.010949141561386,
                       0.010949141561386,
                       0.010949141561386,
                       0.010949141561386};
            return;
        default:
            throw std::runtime_error("Error: Unsupported tetrahedron quadrature polyOrder: " + std::to_string(porder));
    }


}


YAFEL_NAMESPACE_CLOSE